% function matchFFRPath
% 
% Need resolve the model under the zero loadings assumption 

%% Name of the observable that must match 
minforc.nameObs={'FFR','FFR Lead 1','FFR Lead 2','FFR Lead 3','FFR Lead 4','FFR Lead 5',...
    'FFR Lead 6'};  %,'FFR Lead 7','FFR Lead 8','FFR Lead 9','FFR Lead 10'}
%% Name of the shocks that can use 
minforc.nameShocks={'ID Signal 1','ID Signal 2','ID Signal 3','ID Signal 4',...
    'ID Signal 5','ID Signal 6'}; %'Factor 1','Factor 2'
minforc.signalEnd=5; 

%% Number of replications 
minforc.Nrep=10; 

minforc.signalShut=(minforc.signalEnd+1:addsol.NLeadsFFRObs); 
minforc.loadShut  =zeros(length(minforc.signalShut),1); 
minforc.shockShut =zeros(length(minforc.signalShut),1); 

minforc.signalShutPos=[]; 
for ii=1:length(minforc.signalShut); 
    tempstr=strfind(stateNames,['Signal ',num2str(minforc.signalShut(ii))]); 
    tempStru=findNonEmptyCell(tempstr); 
    minforc.signalShutPos=[minforc.signalShutPos;tempStru.indNonEmpty(:)]; 
    minforc.loadShut(ii)=cellposition(['mathcal{B}^{S2}_{',num2str(minforc.signalShut(ii)),'}'],parnames); 
    minforc.shockShut(ii)=cellposition(['ID Signal ',num2str(minforc.signalShut(ii))],shockNames); 
end 
clear temp*; 

parCounter=parMode; 
parCounter(minforc.loadShut)=0; 
disp('states to shut down'); 
stateNames(minforc.signalShutPos)
quer('c'); 


%% Solve model, obtain matrices 
[Gbase, Rbase, Cbase, eu, SDXbase, Zbase]=feval(funcmod,parMode,solveopt,addsol);

[Gcount, Rcount, Ccount, eu, SDXcount, Zcount]...
    =feval(funcmod,parCounter,solveopt,addsol);


minforc.G=Gcount(:,:,end); 
minforc.R=Rcount(:,:,end); 
minforc.Z=Zcount(:,:,end); 
minforc.C=Ccount(:,end); 

%% Positions, Dimensions and data 
minforc.posObs=cellposition(minforc.nameObs,Znames); 
minforc.posInnov=cellposition(minforc.nameShocks,shockNames); 
minforc.data=Y(end,:)'; 
minforc.sTMinusOne=KFStru.smoothSt(end-1,:); 
minforc.innovT=KFStru.innovations(end,:); 
minforc.Nx=size(SDX,1); 
minforc.Nmatch=length(minforc.nameObs); 
minforc.Nfree=length(minforc.nameShocks); 

%% Set desired signals to zero 
if isempty(minforc.signalShutPos)==false; 
    minforc.sTMinusOne(minforc.signalShutPos)=0; 
end 
minforc.innovT(minforc.shockShut)=0; 

%% Initial Guess 
x0=minforc.innovT(minforc.posInnov); 

fval=minObsDiscrep(x0,minforc); 

minproblem.optimset=optimset('MaxIter',1000,'FunValCheck','On','Display','Iter'); 
minproblem.x0=x0; 
minproblem.maxIter=10000; 
minproblem.TolFun=1e-8; 
minproblem.TolX=1e-8; 
outSt=minObslsqnonlin(minproblem,minforc); 

%%
minforc.innovT(minforc.posInnov)=outSt.xstar;
%Generate one step ahead forecast
sT=minforc.G*minforc.sTMinusOne(:)+minforc.R*minforc.innovT(:);
yT=minforc.Z*(minforc.C+sT);
%Compute the differences
fval(:)=minforc.data(minforc.posObs)-yT(minforc.posObs);

Nforc=10; 
tauVec=addsol.tauVec(end)*ones(Nforc,1); 

%% Counterfactual Forecast 
etaMatCount=zeros(Nforc,minforc.Nx); 
etaMatCount(1,:)=(minforc.innovT(:))'; 
[stforcCount,yforcCount]=forecastLoop(Gcount,Rcount,Zcount,Ccount,tauVec,minforc.sTMinusOne(:),...
    etaMatCount,Nforc); 

%% Baseline Forecast 
etaMat=zeros(Nforc,minforc.Nx); 
etaMat(1,:)=KFStru.innovations(end,:)'; 
[stforcBase,yforcBase]=forecastLoop(Gbase,Rbase,Zbase,Cbase,tauVec,KFStru.smoothSt(end-1,:),...
    etaMat,Nforc); 

disp('Checking forecast is correct'); 
toldif=1e-5; 
maxdif=comparemat(yforcBase(1,minforc.posObs),yforcCount(1,minforc.posObs)); 
if maxdif > toldif 
    error('Counterfactual does not recover the obsservables in minforc'); 
end 

%minproblem.maxIter=5000; 
%fstar=lsqnonlin(@(x) minObsDiscrep(x,minforc),x0,minproblem.optimset);
%fstar=lsqnonlin(minproblem,minforc); 